library(tidyverse) #ggplot2, dplyr, tibble, etc.
library(plotly) #use ggplotly(*plotName*) for interactive plots with scoll-over IDs
library(knitr)
library(ggfortify)
library(Rtsne)
library(usethis) #git integration
library(glmmTMB)
library(lme4)
library(ggpubr) #ggarrange() figure wraps
library(plyr)
library(ggeffects) #ggpredict
library(rstanarm) #stan models
library(shinystan) #stan model evaluation >launch_shinystan_demo()
library(loo) #loo() to compare fits between bayesian models
library(MuMIn) #dredge, model averaging
library(RColorBrewer) ; brewer.pal(n=8, name = "Dark2")#; display.brewer.pal(n = 8, name = "Dark2")
## [1] "#1B9E77" "#D95F02" "#7570B3" "#E7298A" "#66A61E" "#E6AB02" "#A6761D"
## [8] "#666666"
theme_set(theme_classic())
df <- data.frame(read.csv("Nestlings1989-2018_no2011_WITH_AdultMeasurements_2001-2015.csv", h=T))
summary(df[,c(2:23)])
## year population species lifeStage
## Min. :1989 Length:4372 Length:4372 Length:4372
## 1st Qu.:1995 Class :character Class :character Class :character
## Median :2002 Mode :character Mode :character Mode :character
## Mean :2002
## 3rd Qu.:2009
## Max. :2018
##
## id age nest nestName
## Length:4372 Min. : 0.00 Length:4372 Length:4372
## Class :character 1st Qu.: 10.00 Class :character Class :character
## Mode :character Median : 24.00 Mode :character Mode :character
## Mean : 28.77
## 3rd Qu.: 27.00
## Max. :22014.00
## NA's :336
## sexSummary tarsus sevenPrimary tailLength
## Length:4372 Min. : 0.00 Min. : 0.0 Min. : 0.00
## Class :character 1st Qu.: 43.80 1st Qu.: 27.0 1st Qu.: 7.00
## Mode :character Median : 57.00 Median : 97.0 Median : 47.00
## Mean : 49.43 Mean : 92.5 Mean : 49.82
## 3rd Qu.: 59.90 3rd Qu.: 121.0 3rd Qu.: 69.00
## Max. :144.00 Max. :14550.0 Max. :195.00
## NA's :30 NA's :100 NA's :123
## calcAge billNT billD billW
## Min. : 0.00 Min. : -0.1224 Min. : 0.00 Min. : 0.00
## 1st Qu.: 13.00 1st Qu.: 16.1000 1st Qu.:11.10 1st Qu.: 11.60
## Median : 24.98 Median : 21.4000 Median :12.40 Median : 12.60
## Mean : 29.91 Mean : 19.8106 Mean :11.95 Mean : 12.14
## 3rd Qu.: 26.31 3rd Qu.: 23.4000 3rd Qu.:13.20 3rd Qu.: 13.40
## Max. :29344.26 Max. :158.0000 Max. :20.40 Max. :112.60
## NA's :551 NA's :48 NA's :60 NA's :60
## TEC head skull wingChord
## Min. : 0.00 Min. : 0.00 Min. :-52.80 Min. : 0.0
## 1st Qu.:24.90 1st Qu.: 60.70 1st Qu.: 35.60 1st Qu.:103.0
## Median :33.20 Median : 73.90 Median : 40.30 Median :157.0
## Mean :30.29 Mean : 67.64 Mean : 37.33 Mean :142.2
## 3rd Qu.:36.00 3rd Qu.: 78.10 3rd Qu.: 42.50 3rd Qu.:181.0
## Max. :58.50 Max. :101.70 Max. : 93.90 Max. :329.0
## NA's :53 NA's :59 NA's :55 NA's :356
## exp7Primary weight
## Min. : 0.00 Min. : 0.0
## 1st Qu.: 39.00 1st Qu.:200.0
## Median : 56.00 Median :345.0
## Mean : 58.26 Mean :291.7
## 3rd Qu.: 73.00 3rd Qu.:400.0
## Max. :2038.00 Max. :615.0
## NA's :2673 NA's :120
Add a new column with the number of occurrences of each individual in the dataset.
df$indivCounts <- as.numeric(ave(df$id, df$id, FUN = length))
#sorted by ascending year and lifeStage (add sorting conditions with: ', variableName')
df <- df %>% arrange(year, lifeStage)
#make lifeStage a factor
df$lifeStage <- as.factor(df$lifeStage)
#check n for both life stages
length(df$lifeStage[df$lifeStage == "nestling"]) #4086 nestlings
## [1] 4086
length(df$lifeStage[df$lifeStage == "adult"]) #286 adults
## [1] 286
#create column with numbered observations for reference
df$obs <- seq(1,length(df$nestlingMeasurementID))
#make sex a factor
df$sexSummary <- as.factor(df$sexSummary)
levels(df$sexSummary)
## [1] "F" "M"
#only individuals with multiple measurements (mm)
mm <- subset(df, df$indivCounts>1)
#only nestlings (ne)
ne <- subset(df, df$lifeStage == "nestling")
#only adults (ad)
ad <- subset(df, df$lifeStage == "adult")
#subset morphometric data from .ne
nestlMorph.ne <- ne[,c(1:23)] #n=4086
#drop NAs
nestlMorph.ne <- drop_na(nestlMorph.ne)#n=1413
nestlMorph.ne$sexSummary <- as.factor(nestlMorph.ne$sexSummary)
class(nestlMorph.ne$sexSummary)
## [1] "factor"
levels(nestlMorph.ne$sexSummary)
## [1] "F" "M"
#filter ages (also reduced by NA removal for unknown sexes)
nestlMorph.ne <- nestlMorph.ne %>% filter(between(age, 23, 27))
#create nestling bill morphometric dataframe
nestlBills.ne <- nestlMorph.ne[,c(15:18)]
tarsusByWeight <- nestlMorph.ne %>%
ggplot(aes(x=weight,y=tarsus,color=sexSummary,label=id))+
geom_point()
ggplotly(tarsusByWeight)
skullByWeight <- nestlMorph.ne %>%
ggplot(aes(x=weight, y=skull, color = sexSummary, label = id))+
geom_point()
ggplotly(skullByWeight)
billPCA <- prcomp(nestlBills.ne)
summary(billPCA)
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 2.5344 0.9077 0.76513 0.5385
## Proportion of Variance 0.7908 0.1014 0.07207 0.0357
## Cumulative Proportion 0.7908 0.8922 0.96430 1.0000
billPCA
## Standard deviations (1, .., p=4):
## [1] 2.5344278 0.9077234 0.7651260 0.5384864
##
## Rotation (n x k) = (4 x 4):
## PC1 PC2 PC3 PC4
## billNT -0.5308950 0.08039039 -0.82350303 -0.1831138
## billD -0.1795232 0.30231519 -0.06245138 0.9340647
## billW -0.1414367 0.90873294 0.24765403 -0.3047419
## TEC -0.8160377 -0.27631026 0.50656694 -0.0335406
billPCA.ne <- tibble(PC1=billPCA$x[,1],
PC2=billPCA$x[,2],
sex=nestlMorph.ne$sexSummary,
id=nestlMorph.ne$id)
billPCA_plot <- billPCA.ne %>%
ggplot(aes(x=PC1,y=PC2,color=sex,label=id))+
geom_point()
ggplotly(billPCA_plot)
morphPCA <- prcomp(nestlMorph.ne[,11:23])
morphPCA.plot <- tibble(PC1=morphPCA$x[,1],
PC2=morphPCA$x[,2],
sex=nestlMorph.ne$sexSummary,
id=nestlMorph.ne$id)
summary(morphPCA)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 53.7309 26.2520 5.45909 4.19985 2.75474 2.28349 1.85734
## Proportion of Variance 0.7922 0.1891 0.00818 0.00484 0.00208 0.00143 0.00095
## Cumulative Proportion 0.7922 0.9813 0.98949 0.99433 0.99642 0.99785 0.99879
## PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 1.64247 0.87822 0.66903 0.50423 0.47736 2.94e-15
## Proportion of Variance 0.00074 0.00021 0.00012 0.00007 0.00006 0.00e+00
## Cumulative Proportion 0.99953 0.99974 0.99987 0.99994 1.00000 1.00e+00
morphPCA
## Standard deviations (1, .., p=13):
## [1] 5.373088e+01 2.625199e+01 5.459091e+00 4.199851e+00 2.754744e+00
## [6] 2.283486e+00 1.857336e+00 1.642469e+00 8.782189e-01 6.690293e-01
## [11] 5.042252e-01 4.773575e-01 2.940383e-15
##
## Rotation (n x k) = (13 x 13):
## PC1 PC2 PC3 PC4 PC5
## tarsus -0.038179895 -0.0041571604 0.005423045 0.120275476 0.366072576
## sevenPrimary -0.180750280 0.5007468830 0.142762776 0.311802527 -0.718103471
## tailLength -0.169665443 0.3659809253 0.555711993 -0.711826105 0.121406957
## calcAge -0.009755084 0.0292277133 0.010712123 -0.013512652 -0.014560321
## billNT -0.018667171 0.0191472815 0.005238650 0.020614391 0.086587176
## billD -0.009043039 0.0009808141 -0.005119253 0.008156555 -0.007620977
## billW -0.005325301 -0.0001540695 0.006722487 0.015456294 0.002797935
## TEC -0.024821963 0.0234086209 -0.003332407 0.012766991 0.107772601
## head -0.045274211 0.0406639533 -0.027252951 0.096631723 0.205755069
## skull -0.020452248 0.0172553324 -0.023920545 0.083864732 0.097982468
## wingChord -0.212187158 0.4938662426 0.197633543 0.511596188 0.506516463
## exp7Primary -0.131362064 0.4958120660 -0.791512479 -0.319731594 0.077758697
## weight -0.933363781 -0.3490539194 -0.060687124 -0.014515091 -0.040588008
## PC6 PC7 PC8 PC9
## tarsus -0.3839684220 0.139908741 0.823167452 -0.0593058494
## sevenPrimary -0.2437940532 0.020134803 0.150395150 0.0096874574
## tailLength -0.0561074066 0.054020382 0.004834509 -0.0073901265
## calcAge 0.0006989063 -0.018584139 0.007881656 0.0025541398
## billNT -0.2503866342 -0.179780817 -0.166879870 0.0002973483
## billD -0.0288488520 -0.066656717 -0.048722349 -0.2476934611
## billW -0.0739067474 -0.060008910 -0.086528860 -0.9562300860
## TEC -0.3451388317 -0.674472290 -0.078289287 0.1003405420
## head -0.6229716734 0.016993897 -0.382415366 0.1014748445
## skull -0.2778328417 0.691466187 -0.304126079 0.0011343025
## wingChord 0.3663868389 -0.054623473 -0.128056173 -0.0002542232
## exp7Primary 0.0318983036 0.008833085 0.025970132 -0.0146684138
## weight 0.0365215872 -0.001528300 -0.006702304 0.0042229533
## PC10 PC11 PC12 PC13
## tarsus -0.0373501140 0.011609262 0.019037700 1.238924e-16
## sevenPrimary -0.0098638877 0.014317600 -0.021248184 1.119139e-16
## tailLength 0.0117819828 0.025549436 -0.009160179 -7.667362e-17
## calcAge 0.0132875073 -0.755419223 0.653675739 6.533814e-16
## billNT -0.9222981215 -0.095524280 -0.093342817 -7.585426e-17
## billD -0.1091091636 0.626340942 0.725642728 6.659476e-16
## billW 0.0837301651 -0.163407019 -0.187191368 2.443571e-16
## TEC 0.2514661116 0.006843275 -0.014155675 -5.773503e-01
## head 0.2554900619 0.009773461 0.015998322 5.773503e-01
## skull 0.0040239503 0.002930186 0.030153998 -5.773503e-01
## wingChord 0.0089653382 0.003504284 -0.003161622 -4.357466e-17
## exp7Primary -0.0059954773 0.005145221 -0.010006872 7.870149e-17
## weight -0.0001802609 -0.005463276 -0.004859557 -2.439204e-17
indMorphPCA <- cbind(morphPCA$x[,1], nestlMorph.ne$id)
morphPCA <- morphPCA.plot %>%
ggplot(aes(x=PC1,y=PC2,color=sex,label=id))+
geom_point()
ggplotly(morphPCA)
nestlMorph.ne %>%
ggplot(aes(x=sexSummary,y=billNT))+
geom_boxplot()
t.test(nestlMorph.ne$billNT[nestlMorph.ne$sexSummary=="M"], nestlMorph.ne$billNT[nestlMorph.ne$sexSummary=="F"])
##
## Welch Two Sample t-test
##
## data: nestlMorph.ne$billNT[nestlMorph.ne$sexSummary == "M"] and nestlMorph.ne$billNT[nestlMorph.ne$sexSummary == "F"]
## t = 3.7865, df = 318.26, p-value = 0.0001826
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.2874067 0.9091245
## sample estimates:
## mean of x mean of y
## 22.09316 21.49490
nestlMorph.ne %>%
ggplot(aes(x=sexSummary,y=billW))+
geom_boxplot()
t.test(nestlMorph.ne$billW[nestlMorph.ne$sexSummary=="M"], nestlMorph.ne$billW[nestlMorph.ne$sexSummary=="F"])
##
## Welch Two Sample t-test
##
## data: nestlMorph.ne$billW[nestlMorph.ne$sexSummary == "M"] and nestlMorph.ne$billW[nestlMorph.ne$sexSummary == "F"]
## t = 5.4804, df = 353.99, p-value = 8.076e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.3285849 0.6964170
## sample estimates:
## mean of x mean of y
## 12.89538 12.38288
nestlMorph.ne %>%
ggplot(aes(x=sexSummary,y=billD))+
geom_boxplot()
t.test(nestlMorph.ne$billD[nestlMorph.ne$sexSummary=="M"], nestlMorph.ne$billD[nestlMorph.ne$sexSummary=="F"])
##
## Welch Two Sample t-test
##
## data: nestlMorph.ne$billD[nestlMorph.ne$sexSummary == "M"] and nestlMorph.ne$billD[nestlMorph.ne$sexSummary == "F"]
## t = 4.3053, df = 332.26, p-value = 2.197e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.1790865 0.4804187
## sample estimates:
## mean of x mean of y
## 12.70203 12.37227
nestlMorph.ne %>%
ggplot(aes(x=sexSummary,y=TEC))+
geom_boxplot()
t.test(nestlMorph.ne$TEC[nestlMorph.ne$sexSummary=="M"], nestlMorph.ne$TEC[nestlMorph.ne$sexSummary=="F"])
##
## Welch Two Sample t-test
##
## data: nestlMorph.ne$TEC[nestlMorph.ne$sexSummary == "M"] and nestlMorph.ne$TEC[nestlMorph.ne$sexSummary == "F"]
## t = 4.9907, df = 313.39, p-value = 9.994e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.6710488 1.5445208
## sample estimates:
## mean of x mean of y
## 34.53278 33.42500
df %>%
ggplot(aes(x=lifeStage,y=tarsus,color=sexSummary))+
geom_boxplot()
## Warning: Removed 30 rows containing non-finite values (stat_boxplot).
df %>%
ggplot(aes(x=lifeStage,y=billW,color=sexSummary))+
geom_boxplot()
## Warning: Removed 60 rows containing non-finite values (stat_boxplot).
df %>%
ggplot(aes(x=lifeStage,y=billD,color=sexSummary))+
geom_boxplot()
## Warning: Removed 60 rows containing non-finite values (stat_boxplot).
df %>%
ggplot(aes(x=lifeStage,y=billNT,color=sexSummary))+
geom_boxplot()
## Warning: Removed 48 rows containing non-finite values (stat_boxplot).
#create morphometric matrix
nestlMorph.mat <- as.matrix(nestlMorph.ne[,-c(1:10,22)])
#run tSNE
nestlMorph.tsne <- Rtsne(nestlMorph.mat)
nestlMorphTsne.plot <- data.frame(x=nestlMorph.tsne$Y[,1],
y=nestlMorph.tsne$Y[,2],
nest=nestlMorph.ne$nest,
population=nestlMorph.ne$population,
sex=nestlMorph.ne$sexSummary,
year=nestlMorph.ne$year)
ggplot(nestlMorphTsne.plot, aes(x=x,y=y,color=as.factor(year)))+
geom_point(size = 2)+
xlab("Dimension 1")+
ylab("Dimension 2")+
theme(#legend.position="none",
legend.title = element_text(size = 18, face = "bold"),
axis.title.x = element_text(size = 18, face = "bold"),
axis.text = element_text(size = 16),
axis.title.y = element_text(size = 18, face = "bold"))
cAge_morphs.ne <- ne %>% filter(between(calcAge, 1,50))
cAge_morphs.ne$year <- as.factor(cAge_morphs.ne$year)
#filter out weights <=1
cAge_morphs.filter.ne <- cAge_morphs.ne %>% filter(weight >=1)
summary(cAge_morphs.ne$calcAge)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 13.00 24.98 20.11 26.28 50.00
summary(cAge_morphs.filter.ne$weight)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.0 240.0 350.0 304.3 397.0 500.0
#ageWtCheck <- as.data.frame(cbind(Wt = cAge_morphs.ne$weight, cAge = cAge_morphs.ne$calcAge))
cAge_morphs.filter.ne %>%
ggplot(aes(x = calcAge, y = weight, color = as.factor(year)))+
geom_point(shape = 1)+
theme(legend.position = "none")
wtMidCluster <- cAge_morphs.filter.ne %>% filter(between(calcAge, 24,30))
wtMidCluster %>%
ggplot(aes(x = calcAge, y = weight, color = as.factor(year)))+
geom_point(shape = 1)+
theme(legend.position = "none")
Decimal points come from the formula for age calculation. Should we be rounding to nearest day?
cAge_morphs.ne %>%
ggplot(aes(x = calcAge, y = tailLength, color = as.factor(year)))+
geom_point(shape = 1)+
theme(legend.position = "none")
## Warning: Removed 90 rows containing missing values (geom_point).
tailMidCluster <- cAge_morphs.ne %>% filter(between(calcAge, 24,30))
tailMidCluster %>%
ggplot(aes(x = calcAge, y = tailLength, color = year))+
geom_point(shape = 1)+
theme(legend.position = "none")
## Warning: Removed 8 rows containing missing values (geom_point).
cAge_morphs.ne %>%
ggplot(aes(x = calcAge, y = sevenPrimary, color = as.factor(year)))+
geom_point(shape = 1)+
theme(legend.position = "none")
## Warning: Removed 74 rows containing missing values (geom_point).
primMidCluster <- cAge_morphs.ne %>% filter(between(calcAge, 24,30))
primMidCluster %>%
ggplot(aes(x = calcAge, y = sevenPrimary, color = year))+
geom_point(shape = 1)+
theme(legend.position = "none")
## Warning: Removed 1 rows containing missing values (geom_point).
cAge_morphs.ne %>%
ggplot(aes(x = calcAge, y = tarsus, color = as.factor(year)))+
geom_point(shape = 1)+
theme(legend.position = "none")
## Warning: Removed 13 rows containing missing values (geom_point).
tarsMidCluster <- cAge_morphs.ne %>% filter(between(calcAge, 24,30))
tarsMidCluster.filter <- tarsMidCluster %>% filter(between(tarsus, 20,100))
tarsMidCluster.filter %>%
ggplot(aes(x = calcAge, y = tarsus, color = year))+
geom_point(shape = 1)+
theme(legend.position = "none")